#####
# Script to generate plots of Kessler region entry time from two input dataframes.
#####

# Function to generate summary statistics and plot of Kessler region entry time histogram
ktime_summary_plot <- function(results_summary){

	ktime_plot <- ggplot(data = results_summary) +
		geom_histogram(aes(x=ktime, y=..density..), binwidth=1) +
		geom_vline(aes(xintercept=ktime[ktime_zero])) +
		theme_bw() +
		labs(x="Kessler region entry time", y="Relative frequency", title="") +
		theme(
			axis.text.x = element_text(size=15),
			axis.text.y = element_text(size=15),
			axis.title.x = element_text(size=15),
			axis.title.y = element_text(size=15)
		)

	return(list(summary=results_summary, plot=ktime_plot))
}

# Function to generate summary statistics and plot of Kessler region entry time (y) against elasticity (x)
ktime_elasticity_plot <- function(results_summary){

	min_ktime <- min(results_summary$ktime)
	max_ktime <- max(results_summary$ktime)
	message(min_ktime, max_ktime)

	results_summary$elasticity <- as.numeric(as.character(results_summary$elasticity))
	results_summary$elasticity <- factor(results_summary$elasticity, levels=sort(unique(results_summary$elasticity)))

	results_summary$betaSD <- as.numeric(as.character(results_summary$betaSD))
	results_summary$betaSD <- factor(results_summary$betaSD, levels=sort(unique(results_summary$betaSD)))

	# Plot ktime (y) against returns growth rate (x)
	ktime_returns_plot <- ggplot(results_summary %>% filter(betaDD==326, betaSD==332)) +
		geom_line(aes(x=returns_growth_rate, y=ktime, color=elasticity), linewidth=2) +
		theme_minimal() +
		labs(x="Annual returns growth rate [%]", y="Kessler region entry time [year]", title="", color="Occupancy\nelasticity") +
		theme(
			axis.text.x = element_text(size=15),
			axis.text.y = element_text(size=15),
			axis.title.x = element_text(size=15),
			axis.title.y = element_text(size=15),
			legend.position = c(0.95, 0.95),
			legend.justification = c("right", "top"),
			legend.background = element_rect(fill="white"),
			legend.box.background = element_rect(fill="white"),
        	legend.text = element_text(size=13),
			legend.title = element_text(size=13)
		) +
		scale_color_grey(start=0.5, end=0) +
		geom_hline(yintercept=2023, linetype="dashed") +
		geom_vline(xintercept=2.5, linetype="dashed") +
		geom_vline(xintercept=0) +
		scale_y_continuous(breaks=c(2023, seq(2100, 2300, by=100)), 
						labels=c("2023", seq(2100, 2300, by=100))) +
		labs(x="Annual returns growth rate [%]",
			y="Kessler region entry time [year]", title="", color="Occupancy\nelasticity") +
   		coord_cartesian(ylim=c(2018, 2300))

	
	# Plot ktime (y) against debris-debris fragmentation parameter (x)
	ktime_betaDD_plot <- ggplot(results_summary %>% filter(returns_growth_rate==3, elasticity==0)) +
		geom_line(aes(x=betaDD, y=ktime, color=betaSD), linewidth=2) +
		theme_minimal() +
		theme(
			axis.text.x = element_text(size=15),
			axis.text.y = element_text(size=15),
			axis.title.x = element_text(size=15),
			axis.title.y = element_text(size=15),
			legend.position = c(0.95, 0.95),
			legend.justification = c("right", "top"),
			legend.background = element_rect(fill="white"),
			legend.box.background = element_rect(fill="white"),
        	legend.text = element_text(size=13),
			legend.title = element_text(size=13)
		) +
		scale_color_grey(start=0.5, end=0) +
		geom_hline(yintercept=2023, linetype="dashed") +
		geom_vline(xintercept=0) +
		scale_y_continuous(breaks=c(2023, seq(2100, 2300, by=100)), 
						labels=c("2023", seq(2100, 2300, by=100))) +
		labs(x="New fragments from debris-debris collisions [fragments/collision]", 
			y="Kessler region entry time [year]", title="", color="New fragments from\nsatellite-debris collisions") +
		coord_cartesian(ylim=c(2018, 2300))



	return(list(summary=results_summary, ktime_returns=ktime_returns_plot, ktime_betaDD=ktime_betaDD_plot))
}

# Load libraries
library(data.table)
library(tidyverse)
library(patchwork)

# Make ktime list
ktime_list <- list()
# Load data
filename_vector <- list.files(path="../data/calibration-data/calibrated-sim-output-data/")[-1]
for(file in seq_along(filename_vector)){
    # Read file into data frame
    ktime_list[[file]] <- read.csv(paste0("../data/calibration-data/calibrated-sim-output-data/",filename_vector[file]))

    # Extract first and third numbers from filename
    numbers <- strsplit(filename_vector[file], split='--')[[1]][2]
    numbers <- strsplit(numbers, split='-')[[1]]
    number1 <- as.numeric(numbers[1])
    number3 <- as.numeric(numbers[3])
	number4 <- as.numeric(sub("\\.csv$", "", numbers[4]))

    # Add Kessler indicator variable
    ktime_list[[file]]$Kessler_ind <- as.numeric(ktime_list[[file]]$Kessler > 0)

    # Add Kessler time variable, computed from Kessler_ind
    ktime_list[[file]]$Kessler_time <- ktime_list[[file]]$year[min(which(ktime_list[[file]]$Kessler_ind >= 1))]

    # Add returns growth rate variable -- first number in filename
    ktime_list[[file]]$returns_growth_rate <- number1

    # Add debris-debris fragmentation variable -- third number in filename
    ktime_list[[file]]$betaDD <- number3

    # Add satellite-debris fragmentation variable -- fourth number in filename
    ktime_list[[file]]$betaSD <- number4
}

# Make the big dataframe
k_time_big <- rbindlist(ktime_list)

# Make the summary dataframe
results_summary <- k_time_big %>%
	group_by(returns_growth_rate, sim, betaDD, betaSD) %>%
	summarise(
		elasticity = elasticity[1],
		kessler = Kessler_ind[n()],
		ktime = min(which(Kessler_ind>=1))+2020,
	) %>%
	mutate(
		ktime_zero = which.min(elasticity^2)
	) %>%
	ungroup() %>%
	arrange(elasticity) %>%
	mutate(
		returns_growth_rate = as.numeric(returns_growth_rate)*100,
		elasticity = as.character(elasticity),
		betaDD = as.numeric(betaDD),
		betaSD = as.character(betaSD)
		)

# Simple fix to alleviate problem: when beta_DD = 0, there is only a single steady state; therefore the Kessler_ind variable is always 1, even though Kessler Syndrome is impossible.
results_summary <- results_summary %>%
	mutate(
		kessler = ifelse(betaDD==0, 0, kessler),
		ktime = ifelse(betaDD==0, Inf, ktime)
		) %>%
	filter(ktime!=Inf)

write.csv(results_summary, "../data/calibration-data/calibrated-sim-output-data/results-summary.csv")

ktime_plots <- ktime_elasticity_plot(results_summary)

ktime_returns <- ktime_plots$ktime_returns

ktime_betaDD <- ktime_plots$ktime_betaDD

ggsave(
	filename=paste0("../images/figure-9--kessler-time-returns-and-betaDD.png"),
	plot= ((ktime_betaDD | ktime_returns) + plot_annotation(tag_levels="A")),
	width=6*(25/9),
	height=6,
	units="in",
	bg="white",
	dpi=300
	)
